#include <math.h>
#include <limits>
#include "reg3dvolumedata.h"
#include "cubes.h"

using namespace PBVP;

const int Reg3DVolumeData::vertinfo[2][2][2] = { { { 0, 1 }, { 4, 5 } }, { { 3,
		2 }, { 7, 6 } } };

Reg3DVolumeData::Reg3DVolumeData(unsigned int dim[3]) {
	m_dim[0] = dim[0];
	m_dim[1] = dim[1];
	m_dim[2] = dim[2];
	init();
}

Reg3DVolumeData::Reg3DVolumeData(unsigned int dx, unsigned int dy,
		unsigned int dz) {
	m_dim[0] = dx;
	m_dim[0] = dy;
	m_dim[0] = dz;
	init();
}

Reg3DVolumeData::Reg3DVolumeData(Reg3DVolumeData *pVolData) {
	pVolData->getDimension(m_dim);
	init();
	pVolData->getOrigin(m_orig);
	pVolData->getSpan(m_span);

	setNumOfVariables(pVolData->numOfVariables());

	for (int i = 0; i < numOfVariables(); i++) {
		m_types[i] = pVolData->m_types[i];
		m_VarNames[i] = pVolData->m_VarNames[i];
		m_Gradients[i] = NULL;
		m_MinValues[i] = pVolData->m_MinValues[i];
		m_MaxValues[i] = pVolData->m_MaxValues[i];
		// copy the data array
		void *data;
		switch (pVolData->m_types[i]) {
		case UCHAR_TYPE:
			data = new unsigned char[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].uchar_ptr, numOfVerts()
					* sizeof(unsigned char));
			m_DataPtrs[i].uchar_ptr = (unsigned char*) data;
			break;
		case USHORT_TYPE:
			data = new unsigned short[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].ushort_ptr, numOfVerts()
					* sizeof(unsigned short));
			m_DataPtrs[i].ushort_ptr = (unsigned short*) data;
			break;
		case FLOAT_TYPE:
			data = new float[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].float_ptr, numOfVerts()
					* sizeof(float));
			m_DataPtrs[i].float_ptr = (float *) data;
			break;
		case INT_TYPE:
			data = new int[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].int_ptr, numOfVerts()
					* sizeof(int));
			m_DataPtrs[i].int_ptr = (int *) data;
			break;
		case DOUBLE_TYPE:
			data = new double[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].double_ptr, numOfVerts()
					* sizeof(double));
			m_DataPtrs[i].double_ptr = (double *) data;
			break;
		}
	}
}

Reg3DVolumeData::~Reg3DVolumeData() {
	clear();
}

void Reg3DVolumeData::init() {
	m_NumOfVar = 0;
	m_xy = m_dim[0] * m_dim[1];
}

void Reg3DVolumeData::clear() {
	for (unsigned int i = 0; i < numOfVariables(); i++) {
		DataType type = m_types[i];
		DataPtr data = m_DataPtrs[i];
		switch (type) {
		case UCHAR_TYPE:
			delete[] data.uchar_ptr;
			break;
		case USHORT_TYPE:
			delete[] data.ushort_ptr;
			break;
		case FLOAT_TYPE:
			delete[] data.float_ptr;
			break;
		case INT_TYPE:
			delete[] data.int_ptr;
			break;
		case DOUBLE_TYPE:
			delete[] data.double_ptr;
		}
		if (m_Gradients[i]) {
			delete[] m_Gradients[i];
		}
	}
}

void Reg3DVolumeData::setNumOfVariables(unsigned int n) {
	// clean up any pre-existing data
	if (numOfVariables() > 0) {
		clear();
	}
	m_NumOfVar = n;
	m_DataPtrs.clear();
	m_DataPtrs.resize(n);
	m_types.clear();
	m_types.resize(n);
	m_VarNames.clear();
	m_VarNames.resize(n);
	m_Gradients.clear();
	m_Gradients.resize(n);
	m_MinValues.resize(n);
	m_MaxValues.resize(n);
	for (int i = 0; i < n; i++) {
		m_types[i] = NO_TYPE;
		m_Gradients[i] = NULL;
		m_MinValues[i] = std::numeric_limits<double>::infinity();
		m_MaxValues[i] = -std::numeric_limits<double>::infinity();
	}
}

void Reg3DVolumeData::setDataArray(void *data, unsigned int var,
		PBVP::DataType type) {
	assert(var < numOfVariables());
	m_types[var] = type;
	switch (type) {
	case UCHAR_TYPE:
		m_DataPtrs[var].uchar_ptr = (unsigned char*) data;
		break;
	case USHORT_TYPE:
		m_DataPtrs[var].ushort_ptr = (unsigned short*) data;
		break;
	case FLOAT_TYPE:
		m_DataPtrs[var].float_ptr = (float *) data;
		break;
	case INT_TYPE:
		m_DataPtrs[var].int_ptr = (int *) data;
		break;
	case DOUBLE_TYPE:
		m_DataPtrs[var].double_ptr = (double *) data;
		break;
	}
}

void Reg3DVolumeData::getCellValues(double vals[8], unsigned int i,
		unsigned int j, unsigned int k, unsigned int var /* = 0 */) const {
	unsigned int v = index2vert(i, j, k);
	vals[0] = getValue(v, var);
	vals[1] = getValue(v + 1, var);
	vals[2] = getValue(v + m_xy + 1, var);
	vals[3] = getValue(v + m_xy, var);
	vals[4] = getValue(v + m_dim[0], var);
	vals[5] = getValue(v + m_dim[0] + 1, var);
	vals[6] = getValue(v + m_xy + m_dim[0] + 1, var);
	vals[7] = getValue(v + m_xy + m_dim[0], var);
}

double Reg3DVolumeData::getValue(float coord[3], unsigned int var) const {
	int idx[3];
	float u[2], v[2], w[2];
	double data[8];
	idx[0] = (int) ((coord[0] - m_orig[0]) / m_span[0]);
	idx[1] = (int) ((coord[1] - m_orig[1]) / m_span[1]);
	idx[2] = (int) ((coord[2] - m_orig[2]) / m_span[2]);

	// clamp to boundary
	bool outside = false;
	for (int i = 0; i < 3; i++) {
		if (idx[i] < 0) {
			idx[i] = 0;
			outside = true;
		} else if (idx[i] >= (m_dim[i] - 1)) {
			idx[i] = m_dim[i] - 1;
			outside = true;
		}
	}
	if (outside) {
		return getValue(idx[0], idx[1], idx[2], var);
	}

	getCellValues(data, idx[0], idx[1], idx[2], var);

	u[1] = (coord[0] - m_orig[0]) / m_span[0] - idx[0];
	u[0] = 1 - u[1];
	v[1] = (coord[1] - m_orig[1]) / m_span[1] - idx[1];
	v[0] = 1 - v[1];
	w[1] = (coord[2] - m_orig[2]) / m_span[2] - idx[2];
	w[0] = 1 - w[1];

	//printf("i, j, k = %d %d %d, u v w = %f %f %f\n", idx[0], idx[1], idx[2], u[1], v[1], w[1]);
	double val = 0;
	for (int k = 0; k < 2; k++) {
		for (int j = 0; j < 2; j++) {
			for (int i = 0; i < 2; i++) {
				val += data[vertinfo[k][j][i]] * u[i] * v[j] * w[k];
			}
		}
	}
	return val;
}

void Reg3DVolumeData::getVertGrad(unsigned int i, unsigned int j,
		unsigned int k, double grad[3], unsigned int var) {
	assert(var < numOfVariables());
	unsigned int v = index2vert(i, j, k);
	if (hasGradient(var)) {
		// read from saved gradients
		grad[0] = m_Gradients[var][3* v ];
		grad[1] = m_Gradients[var][3* v + 1];
		grad[2] = m_Gradients[var][3* v + 2];
	} else {
		/*
		 // using central differences
		 if (i==0) {
		 // use right difference
		 grad[0] = getValue(v+1, var) - getValue(v, var);
		 } else if (i == m_dim[0]-1) {
		 // use left difference
		 grad[0] = getValue(v, var) - getValue(v-1, var);
		 } else {
		 // use central difference
		 grad[0] = (getValue(v+1, var) - getValue(v-1, var)) / 2;
		 }

		 if (j == 0) {
		 grad[1] = getValue(v+m_dim[0], var) - getValue(v, var);
		 } else if (j  ==  m_dim[1]-1) {
		 grad[1] = getValue(v, var) - getValue(v-m_dim[0], var);
		 } else {
		 grad[1] = (getValue(v+m_dim[0], var) - getValue(v-m_dim[0], var)) / 2;
		 }

		 if (k == 0) {
		 grad[2] = getValue(v+m_xy, var) - getValue(v, var);
		 } else if (k == m_dim[2]-1) {
		 grad[2] = getValue(v, var) - getValue(v-m_xy, var);
		 } else {
		 grad[2] = (getValue(v+m_xy, var) - getValue(v-m_xy, var)) / 2;
		 }
		 grad[0] /= m_span[0];
		 grad[1] /= m_span[1];
		 grad[2] /= m_span[2];
		 */
		double v[27];
		unsigned int ix[3], iy[3], iz[3];
		int l, m, n, t;

		ix[0] = (i > 0) ? i - 1 : 0;
		ix[1] = i;
		ix[2] = (i + 1 < m_dim[0]) ? i + 1 : i;
		iy[0] = (j > 0) ? j - 1 : 0;
		iy[1] = j;
		iy[2] = (j + 1 < m_dim[1]) ? j + 1 : j;
		iz[0] = (k > 0) ? k - 1 : 0;
		iz[1] = k;
		iz[2] = (k + 1 < m_dim[2]) ? k + 1 : k;

		t = 0;
		for (n = 0; n < 3; n++) {
			for (m = 0; m < 3; m++) {
				for (l = 0; l < 3; l++) {
					v[t] = getValue(ix[l], iy[m], iz[n], var);
					t++;
				}
			}
		}
		grad[0] = grad[1] = grad[2] = 0;
		for (l = 0; l < 27; l++) {
			grad[0] += x_grad_mask[l] * v[l];
			grad[1] += y_grad_mask[l] * v[l];
			grad[2] += z_grad_mask[l] * v[l];
		}
		grad[0] /= m_span[0];
		grad[1] /= m_span[1];
		grad[2] /= m_span[2];
		// handle boundary conditions
		if (i == 0 || i == m_dim[0] - 1)
			grad[0] *= 2;
		if (j == 0 || j == m_dim[1] - 1)
			grad[1] *= 2;
		if (k == 0 || k == m_dim[2] - 1)
			grad[2] *= 2;
	}
}

void Reg3DVolumeData::getVertGradForward(unsigned int i, unsigned int j,
		unsigned int k, double grad[3], unsigned int var /* = 0 */) {
	if (i < m_dim[0] - 1) {
		grad[0] = (getValue(i + 1, j, k, var) - getValue(i, j, k, var))
				/ m_span[0];
	} else {
		grad[0] = (getValue(i, j, k, var) - getValue(i - 1, j, k, var))
				/ m_span[0]; // x boundary
	}
	if (j < m_dim[1] - 1) {
		grad[1] = (getValue(i, j + 1, k, var) - getValue(i, j, k, var))
				/ m_span[1];
	} else {
		grad[1] = (getValue(i, j, k, var) - getValue(i, j - 1, k, var))
				/ m_span[1]; // y boundary
	}
	if (k < m_dim[2] - 1) {
		grad[2] = (getValue(i, j, k + 1, var) - getValue(i, j, k, var))
				/ m_span[2];
	} else {
		grad[2] = (getValue(i, j, k, var) - getValue(i, j, k - 1, var))
				/ m_span[2]; // z boundary
	}
}

void Reg3DVolumeData::getVertGradBackward(unsigned int i, unsigned int j,
		unsigned int k, double grad[3], unsigned int var /* = 0 */) {
	if (i > 0) {
		grad[0] = (getValue(i, j, k, var) - getValue(i - 1, j, k, var))
				/ m_span[0];
	} else {
		grad[0] = (getValue(i + 1, j, k, var) - getValue(i, j, k, var))
				/ m_span[0];
	}
	if (j > 0) {
		grad[1] = (getValue(i, j, k, var) - getValue(i, j - 1, k, var))
				/ m_span[1];
	} else {
		grad[1] = (getValue(i, j + 1, k, var) - getValue(i, j, k, var))
				/ m_span[1];
	}
	if (k > 0) {
		grad[2] = (getValue(i, j, k, var) - getValue(i, j, k - 1, var))
				/ m_span[2];
	} else {
		grad[2] = (getValue(i, j, k + 1, var) - getValue(i, j, k, var))
				/ m_span[2];
	}
}

void Reg3DVolumeData::getQuadENODerivatives(unsigned int i, unsigned int j,
		unsigned int k, double plus[3], double minus[3], unsigned int var) {
	// TODO
	// Five values needed to evaluate Quadratic ENO derivatives
	// Code of getting seven values is adopted from Xuguo's
	double fx[7], fy[7], fz[7];

	getSevenNeighorValues(i, j, k, fx, 0, var);
	getSevenNeighorValues(i, j, k, fy, 1, var);
	getSevenNeighorValues(i, j, k, fz, 2, var);

	// compute 1st and 2nd derivatives at multiple points
	double d1[6], d2[5];
	deriv2SevenValues(fx, d1, d2, m_span[0]);
	// compute Dx^-
	int n = 2;		// the point to be computed is at index 3
	minus[0] = d1[n];	// first order approximation
	double f2 = (fabs(d2[n-1]) < fabs(d2[n]))? d2[n-1]:d2[n];
	minus[0] += f2*(2*(3-n)-1)*m_span[0];
	// compute Dx^+
	n = 3;
	plus[0] = d1[n];
	f2 = (fabs(d2[n-1]) < fabs(d2[n]))? d2[n-1]:d2[n];
	plus[0] += f2*(2*(3-n)-1)*m_span[0];

	deriv2SevenValues(fy, d1, d2, m_span[1]);
	// compute Dy^-
	n = 2;
	minus[1] = d1[n];
	f2 = (fabs(d2[n-1]) < fabs(d2[n]))? d2[n-1]:d2[n];
	minus[1] += f2*(2*(3-n)-1)*m_span[1];
	// compute Dy^+
	n = 3;
	plus[1] = d1[n];
	f2 = (fabs(d2[n-1]) < fabs(d2[n]))? d2[n-1]:d2[n];
	plus[1] += f2*(2*(3-n)-1)*m_span[1];

	deriv2SevenValues(fz, d1, d2, m_span[2]);
	// compute Dz^-
	n = 2;
	minus[2] = d1[n];
	f2 = (fabs(d2[n-1]) < fabs(d2[n]))? d2[n-1]:d2[n];
	minus[2] += f2*(2*(3-n)-1)*m_span[2];
	// compute Dx^+
	n = 3;
	plus[2] = d1[n];
	f2 = (fabs(d2[n-1]) < fabs(d2[n]))? d2[n-1]:d2[n];
	plus[2] += f2*(2*(3-n)-1)*m_span[2];
}

void Reg3DVolumeData::deriv2SevenValues(double vals[7], double d1[6], double d2[5], float dx)
{
	int i;
	float indx1, indx2, indx3;


	indx1 = 1.0/dx;

	// The first order divided difference
	for (i = 1; i < 5; i++) {
		d1[i] = (vals[i+1] - vals[i]) / dx;
	}

	// The second order divided difference
	for (i = 1; i < 4; i++) {
		d2[i] = (d1[i+1] - d1[i])/(2*dx);
	}
}

double Reg3DVolumeData::getOneNeighborAbsMax(int i, int j, int k, unsigned int var /* = 0 */)
{
	double max = fabs(getValue(i, j, k, var));
	for (int ii = std::max<int>(i-1, 0); ii < std::min<int>(i+2, m_dim[0]); ii++)
	{
		if (fabs(getValue(ii, j, k, var)) > max) max = fabs(getValue(ii, j, k, var));
	}
	for (int jj = std::max<int>(j-1, 0); jj < std::min<int>(j+2, m_dim[1]); jj++)
	{
		if (fabs(getValue(i, jj, k, var)) > max) max = fabs(getValue(i, jj, k, var));
	}
	for (int kk = std::max<int>(k-1, 0); kk < std::min<int>(k+2, m_dim[2]); kk++)
	{
		if (fabs(getValue(i, j, kk, var)) > max) max = fabs(getValue(i, j, kk, var));
	}
	return max;
}

double Reg3DVolumeData::getOneNeighborAbsMin(int i, int j, int k, unsigned int var /* = 0 */)
{
	double min = fabs(getValue(i, j, k, var));
	for (int ii = std::max<int>(i-1, 0); ii < std::min<int>(i+2, m_dim[0]); ii++)
	{
		if (fabs(getValue(ii, j, k, var)) < min) min = fabs(getValue(ii, j, k, var));
	}
	for (int jj = std::max<int>(j-1, 0); jj < std::min<int>(j+2, m_dim[1]); jj++)
	{
		if (fabs(getValue(i, jj, k, var)) < min) min = fabs(getValue(i, jj, k, var));
	}
	for (int kk = std::max<int>(k-1, 0); kk < std::min<int>(k+2, m_dim[2]); kk++)
	{
		if (fabs(getValue(i, j, kk, var)) < min) min = fabs(getValue(i, j, kk, var));
	}
	return min;
}

void Reg3DVolumeData::getSevenNeighorValues(unsigned int i, unsigned int j,
		unsigned int k, double vals[7], unsigned int dir, unsigned int var) {
	switch (dir) {
	case 0: // x direction
		if (i >= 3 && i <= m_dim[0] - 4) {
			for (int l = 0; l < 7; l++) {
				vals[l] = getValue(i - 3 + l, j, k, var);
			}
		} else if (i == 2) {
			for (int l = 1; l < 7; l++) {
				vals[l] = getValue(i - 3 + l, j, k, var);
			}
			vals[0] = 2* vals [1] - vals[2];
		} else if (i == 1) {
			for (int l = 2; l < 7; l++) {
				vals[l] = getValue(i - 3 + l, j, k, var);
			}
			vals[1] = 2* vals [2] - vals[3];
			vals[0] = 2* vals [1] - vals[2];
		} else if (i == 0) {
			for (int l = 3; l < 7; l++) {
				vals[l] = getValue(i - 3 + l, j, k, var);
			}
			vals[2] = 2* vals [3] - vals[4];
			vals[1] = 2* vals [2] - vals[3];
			vals[0] = 2* vals [1] - vals[2];
		} else if (i == m_dim[0] - 3) {
			for (int l = 0; l < 6; l++) {
				vals[l] = getValue(i - 3 + l, j, k, var);
			}
			vals[6] = 2* vals [5] - vals[4];
		} else if (i == m_dim[0] - 2) {
			for (int l = 0; l < 5; l++) {
				vals[l] = getValue(i - 3 + l, j, k, var);
			}
			vals[5] = 2* vals [4] - vals[3];
			vals[6] = 2* vals [5] - vals[4];
		} else if (i == m_dim[0] - 1) {
			for (int l = 0; l < 4; l++) {
				vals[l] = getValue(i - 3 + l, j, k, var);
			}
			vals[4] = 2* vals [3] - vals[2];
			vals[5] = 2* vals [4] - vals[3];
			vals[6] = 2* vals [5] - vals[4];
		}
		break;
	case 1: // y direction
		if (j >= 3 && j <= m_dim[1] - 4) {
			for (int l = 0; l < 7; l++) {
				vals[l] = getValue(i, j - 3 + l, k, var);
			}
		} else if (j == 2) {
			for (int l = 1; l < 7; l++) {
				vals[l] = getValue(i, j - 3 + l, k, var);
			}
			vals[0] = 2* vals [1] - vals[2];
		} else if (j == 1) {
			for (int l = 2; l < 7; l++) {
				vals[l] = getValue(i, j - 3 + l, k, var);
			}
			vals[1] = 2* vals [2] - vals[3];
			vals[0] = 2* vals [1] - vals[2];
		} else if (j == 0) {
			for (int l = 3; l < 7; l++) {
				vals[l] = getValue(i, j - 3 + l, k, var);
			}
			vals[2] = 2* vals [3] - vals[4];
			vals[1] = 2* vals [2] - vals[3];
			vals[0] = 2* vals [1] - vals[2];
		} else if (j == m_dim[1] - 3) {
			for (int l = 0; l < 6; l++) {
				vals[l] = getValue(i, j - 3 + l, k, var);
			}
			vals[6] = 2* vals [5] - vals[4];
		} else if (j == m_dim[1] - 2) {
			for (int l = 0; l < 5; l++) {
				vals[l] = getValue(i, j - 3 + l, k, var);
			}
			vals[5] = 2* vals [4] - vals[3];
			vals[6] = 2* vals [5] - vals[4];
		} else if (j == m_dim[1] - 1) {
			for (int l = 0; l < 4; l++) {
				vals[l] = getValue(i, j - 3 + l, k, var);
			}
			vals[4] = 2* vals [3] - vals[2];
			vals[5] = 2* vals [4] - vals[3];
			vals[6] = 2* vals [5] - vals[4];
		}
		break;
	case 2: // z direction
		if (k >= 3 && k <= m_dim[2] - 4) {
			for (int l = 0; l < 7; l++) {
				vals[l] = getValue(i, j, k - 3 + l, var);
			}
		} else if (k == 2) {
			for (int l = 1; l < 7; l++) {
				vals[l] = getValue(i, j, k - 3 + l, var);
			}
			vals[0] = 2* vals [1] - vals[2];
		} else if (k == 1) {
			for (int l = 2; l < 7; l++) {
				vals[l] = getValue(i, j, k - 3 + l, var);
			}
			vals[1] = 2* vals [2] - vals[3];
			vals[0] = 2* vals [1] - vals[2];
		} else if (k == 0) {
			for (int l = 3; l < 7; l++) {
				vals[l] = getValue(i, j, k - 3 + l, var);
			}
			vals[2] = 2* vals [3] - vals[4];
			vals[1] = 2* vals [2] - vals[3];
			vals[0] = 2* vals [1] - vals[2];
		} else if (k == m_dim[2] - 3) {
			for (int l = 0; l < 6; l++) {
				vals[l] = getValue(i, j, k - 3 + l, var);
			}
			vals[6] = 2* vals [5] - vals[4];
		} else if (k == m_dim[2] - 2) {
			for (int l = 0; l < 5; l++) {
				vals[l] = getValue(i, j, k - 3 + l, var);
			}
			vals[5] = 2* vals [4] - vals[3];
			vals[6] = 2* vals [5] - vals[4];
		} else if (k == m_dim[2] - 1) {
			for (int l = 0; l < 4; l++) {
				vals[l] = getValue(i, j, k - 3 + l, var);
			}
			vals[4] = 2* vals [3] - vals[2];
			vals[5] = 2* vals [4] - vals[3];
			vals[6] = 2* vals [5] - vals[4];
		}
		break;
	default:
		// should not be here
		assert(0);
	}
}
void Reg3DVolumeData::getVertParDeriv2(unsigned int i, unsigned int j,
		unsigned int k, double par2[6], unsigned int var /* = 0 */) {
	double v[3][3][3];
	unsigned int ix[3], iy[3], iz[3];
	int l, m, n, t;

	ix[0] = (i > 0) ? i - 1 : 0;
	ix[1] = i;
	ix[2] = (i + 1 < m_dim[0]) ? i + 1 : i;
	iy[0] = (j > 0) ? j - 1 : 0;
	iy[1] = j;
	iy[2] = (j + 1 < m_dim[1]) ? j + 1 : j;
	iz[0] = (k > 0) ? k - 1 : 0;
	iz[1] = k;
	iz[2] = (k + 1 < m_dim[2]) ? k + 1 : k;

	t = 0;
	for (n = 0; n < 3; n++) {
		for (m = 0; m < 3; m++) {
			for (l = 0; l < 3; l++) {
				v[n][m][l] = getValue(ix[l], iy[m], iz[n], var);
				t++;
			}
		}
	}

	par2[0] = (v[1][1][2] - 2* v [1][1][1] + v[1][1][0]) / (m_span[0]
			* m_span[0]); // Fxx
	par2[1] = (v[1][2][2] + v[1][0][0] - v[1][0][2] - v[1][2][0]) / (4*
			m_span [0] * m_span[1]); // Fxy
	par2[2] = (v[2][1][2] + v[0][1][0] - v[2][1][0] - v[0][1][2]) / (4*
			m_span [0] * m_span[2]); // Fxz
	par2[3] = (v[1][2][1] - 2* v [1][1][1] + v[1][0][1]) / (m_span[1]
			* m_span[1]); // Fyy
	par2[4] = (v[2][2][1] + v[0][0][1] - v[2][0][1] - v[0][2][1]) / (4*
			m_span [1] * m_span[2]); // Fyz
	par2[5] = (v[2][1][1] - 2* v [1][1][1] + v[0][1][1]) / (m_span[2]
			* m_span[2]); // Fzz
}

void Reg3DVolumeData::getCellGrads(unsigned int i, unsigned int j,
		unsigned int k, double grads[8][3], unsigned int var) {
	getVertGrad(i, j, k, grads[0], var);
	getVertGrad(i + 1, j, k, grads[1], var);
	getVertGrad(i + 1, j, k + 1, grads[2], var);
	getVertGrad(i, j, k + 1, grads[3], var);
	getVertGrad(i, j + 1, k, grads[4], var);
	getVertGrad(i + 1, j + 1, k, grads[5], var);
	getVertGrad(i + 1, j + 1, k + 1, grads[6], var);
	getVertGrad(i, j + 1, k + 1, grads[7], var);
}

void Reg3DVolumeData::getVertGrad(float x, float y, float z, double grad[],
		unsigned int var) {
	grad[0] = grad[1] = grad[2] = 0;
	int idx[3];
	float u[2], v[2], w[2];
	double cellGrads[8][3];
	idx[0] = (int) ((x - m_orig[0]) / m_span[0]);
	idx[1] = (int) ((y - m_orig[1]) / m_span[1]);
	idx[2] = (int) ((z - m_orig[2]) / m_span[2]);

	if (idx[0] < 0 || idx[0] >= (m_dim[0] - 1))
		return;
	if (idx[1] < 0 || idx[1] >= (m_dim[1] - 1))
		return;
	if (idx[2] < 0 || idx[2] >= (m_dim[2] - 1))
		return;

	getCellGrads(idx[0], idx[1], idx[2], cellGrads, var);

	u[1] = (x - m_orig[0]) / m_span[0] - idx[0];
	u[0] = 1 - u[1];
	v[1] = (y - m_orig[1]) / m_span[1] - idx[1];
	v[0] = 1 - v[1];
	w[1] = (z - m_orig[2]) / m_span[2] - idx[2];
	w[0] = 1 - w[1];
	for (int k = 0; k < 2; k++) {
		for (int j = 0; j < 2; j++) {
			for (int i = 0; i < 2; i++) {
				for (int n = 0; n < 3; n++) {
					grad[n] += cellGrads[vertinfo[k][j][i]][n] * u[i] * v[j]
							* w[k];
				}
			}
		}
	}
}

void Reg3DVolumeData::calcGradient() {
	for (unsigned int var = 0; var < numOfVariables(); var++) {
		if (hasGradient(var))
			continue;

		// Use Bspline approximation to calc gradients
		float v[27];
		unsigned int ix[3], iy[3], iz[3];
		int i, j, k, l, m, n, t;

		double *p_grad = new double[3* m_dim [0] * m_dim[1] * m_dim[2]];
		m_Gradients[var] = p_grad;
		memset(p_grad, 0, sizeof(double) * 3* m_dim [0] * m_dim[1] * m_dim[2]);
		for (k = 0; k < m_dim[2]; k++)
			for (j = 0; j < m_dim[1]; j++)
				for (i = 0; i < m_dim[0]; i++) {
					ix[0] = (i - 1 >= 0) ? i - 1 : 0;
					ix[1] = i;
					ix[2] = (i + 1 < m_dim[0]) ? i + 1 : i;
					iy[0] = (j - 1 >= 0) ? j - 1 : 0;
					iy[1] = j;
					iy[2] = (j + 1 < m_dim[1]) ? j + 1 : j;
					iz[0] = (k - 1 >= 0) ? k - 1 : 0;
					iz[1] = k;
					iz[2] = (k + 1 < m_dim[2]) ? k + 1 : k;

					t = 0;
					for (n = 0; n < 3; n++) {
						for (m = 0; m < 3; m++) {
							for (l = 0; l < 3; l++) {
								v[t] = getValue(ix[l], iy[m], iz[n], var);
								t++;
							}
						}
					}

					double *ptr = p_grad + 3* index2vert (i, j, k);
					for (l = 0; l < 27; l++) {
						*ptr += x_grad_mask[l] * v[l];
						*(ptr + 1) += y_grad_mask[l] * v[l];
						*(ptr + 2) += z_grad_mask[l] * v[l];
					}
					*ptr /= m_span[0];
					*(ptr + 1) /= m_span[1];
					*(ptr + 2) /= m_span[2];

					// handle boundary conditions
					if (i == 0 || i == m_dim[0] - 1)
						*ptr *= 2;
					if (j == 0 || j == m_dim[1] - 1)
						*(ptr + 1) *= 2;
					if (k == 0 || k == m_dim[2] - 1)
						*(ptr + 2) *= 2;
				}
	}
}

double Reg3DVolumeData::getMinValue(unsigned int var) {
	if (m_MinValues[var] == std::numeric_limits<double>::infinity()) {
		calcMinMax();
	}
	return m_MinValues[var];
}

double Reg3DVolumeData::getMaxValue(unsigned int var) {
	if (m_MaxValues[var] == -std::numeric_limits<double>::infinity()) {
		calcMinMax();
	}
	return m_MaxValues[var];
}

void Reg3DVolumeData::calcMinMax() {
	for (unsigned int var = 0; var < numOfVariables(); var++) {
		m_MinValues[var] = m_MaxValues[var] = getValue(0, var);
		for (int i = 0; i < numOfVerts(); i++) {
			if (m_MinValues[var] > getValue(i, var))
				m_MinValues[var] = getValue(i, var);
			if (m_MaxValues[var] < getValue(i, var))
				m_MaxValues[var] = getValue(i, var);
		}
	}
}

void Reg3DVolumeData::resampleUChar(unsigned char ucArray[],
		unsigned int dim[3], unsigned int var) {
	unsigned int i, j, k;
	float minExt[3], maxExt[3];
	float len[3]; // length of new span
	float pos[3];
	double minVal, maxVal;

	getMinMaxExt(minExt, maxExt);
	len[0] = (maxExt[0] - minExt[0]) / (dim[0] - 1);
	len[1] = (maxExt[1] - minExt[1]) / (dim[1] - 1);
	len[2] = (maxExt[2] - minExt[2]) / (dim[2] - 1);
	minVal = getMinValue(var);
	maxVal = getMaxValue(var);
	unsigned int idx = 0;
	for (k = 0; k < dim[2]; k++) {
		for (j = 0; j < dim[1]; j++) {
			for (i = 0; i < dim[0]; i++) {
				pos[0] = m_orig[0] + i * len[0];
				pos[1] = m_orig[1] + j * len[1];
				pos[2] = m_orig[2] + k * len[2];
				double val = getValue(pos, var);
				unsigned char vc;
				if (m_types[var] == UCHAR_TYPE) {
					vc = (unsigned char) val;
				} else {
					vc = (unsigned char) ((val - minVal) * 255 / (maxVal
							- minVal));
				}
				// unsigned int idx = i + j*dim[0] + k*dim[0]*dim[1];
				ucArray[idx] = vc;
				idx++;
			}
		}
	}
}

void Reg3DVolumeData::resampleUCharSimple(unsigned char ucArray[],
		unsigned int dim[], unsigned int var) {
	assert(dim[0] >= width() && dim[1] >= height() && dim[2] >= depth());
	double minVal = getMinValue(var);
	double maxVal = getMaxValue(var);
	unsigned int i, j, k;

	for (k = 0; k < depth(); k++) {
		for (j = 0; j < height(); j++) {
			for (i = 0; i < width(); i++) {
				double val = getValue(i, j, k, var);
				unsigned char vc = (unsigned char) ((val - minVal) * 255
						/ (maxVal - minVal));
				unsigned int idx = i + j * dim[0] + k * dim[0] * dim[1];
				ucArray[idx] = vc;
			}
		}
	}
}

void Reg3DVolumeData::mergeData(PBVP::Reg3DVolumeData *pVolData) {
	unsigned int newSize = m_NumOfVar + pVolData->numOfVariables();
	m_types.resize(newSize);
	m_DataPtrs.resize(newSize);
	m_VarNames.resize(newSize);
	m_Gradients.resize(newSize);
	m_MinValues.resize(newSize);
	m_MaxValues.resize(newSize);

	for (int i = 0; i < pVolData->numOfVariables(); i++) {
		m_types[m_NumOfVar + i] = pVolData->m_types[i];
		m_VarNames[m_NumOfVar + i] = pVolData->m_VarNames[i];
		m_Gradients[m_NumOfVar + i] = NULL;
		m_MinValues[m_NumOfVar + i] = pVolData->m_MinValues[i];
		m_MaxValues[m_NumOfVar + i] = pVolData->m_MaxValues[i];
		// copy the data array
		void *data;
		switch (pVolData->m_types[i]) {
		case UCHAR_TYPE:
			data = new unsigned char[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].uchar_ptr, numOfVerts()
					* sizeof(unsigned char));
			m_DataPtrs[m_NumOfVar + i].uchar_ptr = (unsigned char*) data;
			break;
		case USHORT_TYPE:
			data = new unsigned short[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].ushort_ptr, numOfVerts()
					* sizeof(unsigned short));
			m_DataPtrs[m_NumOfVar + i].ushort_ptr = (unsigned short*) data;
			break;
		case FLOAT_TYPE:
			data = new float[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].float_ptr, numOfVerts()
					* sizeof(float));
			m_DataPtrs[m_NumOfVar + i].float_ptr = (float *) data;
			break;
		case INT_TYPE:
			data = new int[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].int_ptr, numOfVerts()
					* sizeof(int));
			m_DataPtrs[m_NumOfVar + i].int_ptr = (int *) data;
			break;
		case DOUBLE_TYPE:
			data = new double[numOfVerts()];
			memcpy(data, pVolData->m_DataPtrs[i].double_ptr, numOfVerts()
					* sizeof(double));
			m_DataPtrs[m_NumOfVar + i].double_ptr = (double *) data;
			break;
		}
	}
	m_NumOfVar = newSize;
}
